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The continuum theory of partially fluidized shear granular flows is tested and calibrated using two 
f"^ | dimensional soft particle molecular dynamics simulations. The theory is based on the relaxational 

£\j , dynamics of the order parameter that describes the transition between static and flowing regimes of 

granular material. We define the order parameter as a fraction of static contacts among all contacts 
between particles. We also propose and verify by direct simulations the constitutive relation based 
on the splitting of the shear stress tensor into a "fluid part" proportional to the strain rate tensor, 
and a remaining "solid part" . The ratio of these two parts is a function of the order parameter. The 
rheology of the fluid component agrees well with the kinetic theory of granular fluids even in the 
dense regime. Based on the hysteretic bifurcation diagram for a thin shear granular layer obtained 
in simulations, we construct the "free energy" for the order parameter. The theory calibrated using 
numerical experiments with the thin granular layer is applied to the surface-driven stationary two 
O , dimensional granular flows in a thick granular layer under gravity. 
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PACS numbers: 46.55.+d, 45.70.Cc, 46.25.-y 



I. INTRODUCTION 



In the last few year there have been many experimental U E H II IS IS and theoretical H H E IH 13 ES 13 

studies that explored a broad range of granular flow conditions from rapid dilute flows to slow dense flows, as well 
as the details of the shear-driven fluidization transition. While dilute granular flows can be well described by the 
kinetic theory of dissipative granular gases [Til ], dense granular flows still present significant difficulty in formulation 
• of a continuous theory. In Ref. |lq ] , Savage proposed a continuum theory for slow dense granular flows based on the 
so-called associated flow rule that relates the strain rate and the shear stress in plastic frictional systems. Averaging 
strain rate fluctuations yields a Bingham-like constitutive relation in which the shear stress has a viscous as well as 
a strain-rate independent parts. According to this theory, the stress and strain rate tensors are always co-axial and, 
furthermore, it also postulates that the viscosity diverges as the density approaches close packing limit. Losert et al. 
(see also 0) proposed a similar hydrodynamic model based on a Newtonian stress-strain constitutive relation 
with density dependent viscosity without strain-rate independent component. As observed in Ref. Q, the ratio of 
C^l 1 the full shear stress to the strain rate diverges at the fluidization threshold. This was also interpreted in Ref. Q 
as a divergence of the viscosity coefficient when the volume fraction approaches the randomly packed limit. This 
description works only in a fluidized state and can not properly account for hysteretic phenomena in which static and 
fluidized states co-exist under the same external load, such as stick-slip oscillations Q, avalanching or shear band 
formation. 

In many granular flows of interest static and dynamics regions co-exist under the same external load conditions. 
Examples of such hysteretic phenomena include stick-slip oscillations 0], avalanching or shear band formation. 
This calls for a unified theory which would be applicable both in the flowing regime and in the static regime. 

In our recent papers 0, U| we proposed a different approach based on the order parameter description of the 
O ■ granular matter. The value of the order parameter specifies the ratio between static and fluid parts of the stress 
tensor. The order parameter was assumed to obey dissipative dynamics governed by a free energy functional with 
. two local minima. This description based on the separation of static and fluid components of the shear stress, calls 
for an alternative definition of viscosity as a ratio of the fluid part of the shear stress to the strain rate. Since the 
fluid shear stress vanishes together with the strain rate, the viscosity coefficient in our theory is expected to remain 
finite at the fluidization threshold. We assumed the simplest Newtonian friction law, so the viscosity coefficient is a 
constant. This model yielded a good qualitative description of many phenomena occurring in granular flows, such as 
hysteretic transition to chute flow, stick-slip regime of a driven near-surface flow, structure of avalanches in shallow 
chute flows, etc. 

However, several important issues have not been addressed: mesoscopic definition of the order parameter, quantita- 
tive specification of the order parameter dynamics and the constitutive relation. In this paper we set out to perform 
two dimensional (2D) molecular dynamics simulations which should provide us with the way to achieve these goals. 
We classify all contacts as either "fluid-like" or "solid-like" and define the order parameter as a mesoscopic space- 
time average fraction of solid-like contacts. Using this order parameter we obtain the constitutive relation and the 
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relaxational dynamics of the order parameter directly from simulations of granular flow in a thin Couette geometry 
at zero gravity. Preliminary account of our results is presented in Rcf . 20] . 

The paper is organized as follows. Section ^ outlines the standard granular hydrodynamics theory based on the 
kinetic theory of dissipative granular gases. Section IlIII introduces the continuum description of partially fluidized 
flows based on the relaxational dynamics of the order parameter. In Section IIVI we describe our 2D molecular 
dynamics simulations and define measurement protocol for the order parameter and fluid and solid components of 
the stress tensor. In Section we study Couette flow in a thin granular layer to obtain the free energy controlling 
the relaxational order parameter dynamics and to extract the constitutive relations. In Section Ivll the obtained set of 
equations is used to calculate the stress and velocity distributions in a different system, a thick granular layer under 
non-zero gravity driven by a moving heavy upper plate. 



II. GRANULAR HYDRODYNAMICS 



In this Section we outline the standard continuum description of granular flows based on continuity equations for 
mass, momentum and fluctuation kinetic energy (or "granular temperature"). This description is usually applied 
to dilute granular gases where it can be rigorously derived from the kinetic theory [l5| . although slightly modified 
hydrodynamics based on kinetic theory often works reasonably well for relatively dense flows, even though the kinetic 
theory itself is not applicable to these conditions. 

The mass, momentum and energy conservation equations have the usual form 



Dv 

~Dt 
Du 

' Dt 
DT 

' ~Dt 



-vV ■ u, 

-V • <T + ^g, 

er : 7 — V • q — e, 



(1) 
(2) 
(3) 



where v is the density, u is the velocity field, T = (uiiii)/2 is the granular temperature, D/Dt = dt + (u • V) is the 
material derivative, g is the gravity acceleration, rj a p is the stress tensor, q is the energy flux vector, 7 Q ^ = d a ui3+dpu a 
is the strain rate tensor, and e is the energy dissipation rate. 

These three equations have to be supplemented by the constitutive relations for the stress tensor er, energy flux q, 
and the energy dissipation rate e. For dilute systems, a linear relations between stress er and strain rate 7 is obtained, 



cr a f3 = [p + (m - X)Tij]S af 3 - fi%p, 
q = -kVT, 

In the kinetic theory of granular gases 0, these equations are closed with the following equation of state 

AvT 

P = -a l + (l + e)G(i/ , 
■Ear 



and the expressions for the shear and bulk viscosities 
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2TT 1 / 2 dG(ls) 



1 + 2G(v) +(!+-) G{v) 



the thermal conductivity 



7r 1 /2dG(^) 



1 *Gv>)~\\ + i)Gwr 



(4) 
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(6) 

(7) 
(8) 

(9) 



and the energy dissipation rate 



(10) 
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Here < e < 1 is restitution coefficient and d is particle diameter. The function G(v) which enters these relations is 
the spatial particle-particle correlation function, and for a dilute 2D gas of elastic hard disks was derived by Carnahan 
and Starling jlj, 

r (,a v{1 7t//16) fin 

This formula is expected to work for densities roughly below 0.7. For high density granular gases, this function has 
been calculated using free volume theory p^. 

GFV = (1 + e) [(uImW-1] (12) 

where v c « 0.82 is the density of the random close packing limit. Luding [24[ proposed a global fit 

G L = G CS + (1 + exp(-(iy - i^/mo))- 1 )^^ - G CS ) (13) 

with empirically fitted parameters vq 0.7 and mo ~ 10 -2 . However, even with this extension, the continuum theory 
comprised of Eas.l(T))- (|10|l cannot describe the force chains which transmit stress via persistent contacts remaining in 
the dense granular flows, as well as the transition from solid to static regimes and coexisting solid and fluid phases. 



III. ORDER PARAMETER DESCRIPTION OF PARTIALLY FLUIDIZED GRANULAR FLOWS 

In this Section we review briefly our continuum theory [Tsl Il9j which provides an alternative approach to the 
formulation of the constitutive relations in partially fluidized granular flows. In the dense flow regime, the granular 
matter can be considered incompressible, so Eq. JIJ can be replaced by V • u = and the density v = v c . This also 
allows us to drop the energy equation © and the equation of state ©, as for v — > u c , G(v) — > oo and T — > 0, so that 
G{v)T — > const. This of course leads to the familiar divergence of the viscosity coefficient n oc G(i/)T 1 / 2 . 

Next, we separate the stress tensor er into two parts, a static (contact) part <x s , and a fluid part cr? . The latter is 
assumed to take a purely Newtonian form 

a Lf3 = P/<W ~ P>fiaf) (14) 

where pf is the "partial" fluid pressure, ///is the viscosity coefficient associated with the fluid stress tensor which is 
different from //o introduced for the full stress tensor. As we shall see in the following, unlike //q, /// does not diverge 
as v — > v c . In our original model |l8lll9| we simply assumed //,/ = const. 

We postulated [2!| that the fluid part of the off-diagonal components of the stress tensor is proportional to the 
off-diagonal components of the full stress tensor with the proportionality coefficient being a function of the order 
parameter p, 

°{x = q{p)°yx- (15) 

This assumption stipulates the equation for the static part of the off-diagonal stress components, 

o- s yx = (l-q(p))a yx . (16) 

Both fluid and solid parts of the stress tensor are assumed symmetric, = a['y . This assumption is confirmed by our 
numerical simulations (see below). We choose a fixed range for the order parameter such that it is zero in a completely 
fluidized state and one in a completely static regime. Thus the function q{p) has the property q(0) — 1, q(l) = 0. In 
Refs. p!3.ll9| for simplicity we took q(p) = 1 — p. 

A similar relationship can be postulated for the diagonal terms of the stress tensor, 

°"4 = Qx(p)o- xx , <J f yy = q y {p)o-yy, (17) 
°xx = (1 - qx(p))o- xx , a s vv = (1 - q v (p))o- yy , (18) 

where the scaling functions q x ^ y (p) can differ from q(p). 

Combining Eas. (|14|) - (|18[l . we obtain the constitutive relation in the closed form, 



Ca/3 =p f 5 a 0/q a (p) - p f j al3 /q(p) 



(19) 
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where a,/3 £ {x, y}. 

The order parameter itself was not related to any microscopic properties of granular assemblies in Refs. [l8L fli^. 
We simply assumed that because of strong dissipation in dense granular flows it has purely relaxational dynamics 
controlled by the Ginzburg-Landau equation, 

(20, 

Here D is the diffusion coefficient and J-(p) is the free energy density which was assumed to have a quartic polynomial 
form to account for the bistability near the solid-fluid transition, 

HP)= f P (p-l)(p-S)dp. (21) 



The control parameter S is determined by the stress tensor which in Ref. [la . [19| was taken to be a linear function 
of <f) = max \cr mn /a nn \, where the maximum is sought over all possible orthogonal directions m and n. It is easy to 
see that in the interval < S < 1 Eq. (|2U|) has two stable uniform solutions p = 0, 1 corresponding to fluid and solid 
states and one unstable solution p = 8. 

Momentum conservation equation together with Eqs. H19[) - <|21|) represent a closed set of continuous equations 
which after being augmented by appropriate boundary conditions, can describe a variety of interesting granular flows 
such as avalanches in thin chute flows, drum flows, stick-slip oscillation in surface-driven flows, etc. 0,^1- Since 
a rigorous derivation of the continuum model of dense partially fluidized flows from first principles does not seem 
feasible, the assumptions made ad hoc in the formulation of the model l(2*)l. (|19[> - (|21[l have to be checked against 
available experimental and/or numerical data. Unfortunately, it appears to be extremely difficult to directly measure 
the order parameter in the bulk of granular material, as it is determined by subtle changes in the contact fabric. In 
this paper we attempt to extract the properties of the order parameter from molecular dynamics simulations and on 
this basis make a quantitative fit of the order parameter model. 



IV. MOLECULAR DYNAMICS SIMULATIONS 



To model the interaction of individual grains we use the so-called soft-contact approach. The grains are assumed 
to be non-cohesive, dry, inelastic disk-like particles. Two grains interact via normal and shear forces whenever they 
overlap. For the normal impact we employ spring- dashpot model |9j . This model accounts for repulsion and dissipation; 
the repulsive component is proportional to the degree of the overlap, and the velocity dependent damping component 
simulates the dissipation. The model for shear force is based upon the technique developed by Cundall and Strack 
[2^| . It incorporates tangential elasticity and Coulomb laws of friction. The elastic restoring force is proportional to 
the integrated tangential displacement during the contact and limited by the product of the friction coefficient and the 
instantaneous normal force. The grains possess two translational and one rotational degrees of freedom. The motion 
of a grain is obtained by integrating the Newton's equations with the forces and torques produced by its interactions 
with all the neighboring grains and walls of the container. 

Consider a grain i of radius R 1 located at r l moving with translational velocity v 1 and angular velocity a/. This 
grain is in contact with grain j whenever overlap 5 n — R l + R J — \r l — r 4 > 0. The relative velocity at the contact 
point and its normal and tangent components are given by 

v ij = v { -v j + (R i u i +R j u j )t ij , (22) 
v n = (v 1 - v j )n ij , v t = (v l - v j )n ij + {R 1 uj 1 +R j u j ), (23) 

where n IJ = (r l — T l )/r % ^ = (n x ,n y ) is the inward normal to the surface of i at the contact point with j, and the 
direction of the tangent t y = (n y , —n x ) is chosen so that t 4 - 7 x n 1 -? is co-linear with the angular velocity. Then the 
grain i is subjected to the contact force due to the interaction with j with normal and tangential components 

F„ = k n S n - 2~f n m e v n , (24) 

Ft = -sign(5 t ) min(\k t S t \,\fxF n \), 5 t (t) = f v t {r)dr, (25) 

where fc„ jt are the corresponding spring constants, 7„ is the normal damping coefficient, m e = m l m 3 /(rrf + m J ) is 
the reduced mass, p is coefficient of friction, and St is tangential displacement since the moment to of the initial 
contact between i and j. When the static yield criterion Eq. (|25|l is satisfied, the magnitude of S t is adjusted to an 
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instantaneous equilibrium value providing F t = (J,F n . According to the analytical solution of the linear spring-dashpot 
model the coefficient of restitution and the duration of heads-on collision are e = exp(— 7„i n ), t n = n/(k n /m e — 7 2 ) 1 / 2 - 

The advantages and limitations of the employed contact force model were thoroughly studied by a number of 
authors 0, 0, 0] . In fact this is the simplest model which allows us to account for both static and dynamic friction. 

When all forces acting on grain i from other grains, boundaries and perhaps external fields are computed, the 
problem is reduced to the integration of the Newton's equations for translational and rotational degrees of freedom, 

C 

R^Ft, (27) 

c 

■ 2 

where the mass of particle i is denoted by vn l and its moment of inertia is P — l/2m t R l , m l g stands for an external 
gravity field, and the sums in Eqs. i|2ti[l . (|27[l run over a ^ contacts of particle i. 

The results of the simulations reported here are presented in dimensionless form. All quantities are normalized by 
an appropriate combination of the average particle diameter <i, mass m, and gravity g. 

Eqs. P^ )l .p7 )l were integrated using the fifth-order predictor-corrector [23 with a constant time step St. The 
spring constant k n and damping coefficient 7„ were chosen to provide the desired value of the restitution coefficient 
e and guarantee an accurate resolution of an individual collision. Typically, we used e = 0.92, St = 1(T 4 , t c = 505t, 
k n = 2.0 10 5 , 7 „ = 16.7, k s = l/3fc„. 

The computational domain spans L x x L y area, and is periodic in horizontal direction x. Unless otherwise mentioned 
the granulate is slightly polydisperse to avoid crystallization effects @- We assume that the grain diameters are 
uniformly distributed around mean with relative width A r . To provide a link between micro- mechanical quantities 
obtained through simulations and continuous fields, we define the following coarse-graining procedure. Since all 
experiments described below deal with steady quantities, the procedure consists of two steps: space and then time 
averaging. The space averaging of a field x is performed over horizontal bins (along the flow direction) of size 
Vi, f» L x x d and is denoted as (x) in the following. Contributions of those particles which only partially belong to a 
certain bin are weighted by the fraction of their area. After a simulation has reached a steady state, instantaneous 
profiles are averaged over a suitable number of time-snapshots. We shall denote the time-averaging with x. 

For example, steady solid fraction profile is given by 

Kt/) = (K^)), (v( y ,t)) = v b - 1 J2 wi (y) vi > ( 28 ) 

iev b 

where the summation runs over grains at least partially in a bin V& centered at y, V" 1 is the area of grain i, and w l {y) 
is a corresponding fraction of a grain's area within Vb- The coarse-grained velocity field is 

u( y ) ="k^*)>, <u(y,*)> = K^rVM) = Kv.tr'n -1 E w'ww. (29) 

For simplicity we extend the space-time averaging technique described above for other quantities, such as the stress 
tensor 

a af3 (y) = (a aP (y,t)), (a a p(y,t)) = X>M^ + < m ^4>> ( 30 ) 

where a,/3= {x,y}, = r lc ■ e a , FJf — F zc ep, v z a — v l a — {v l a }. The stress tensor in Eq. lp!H| has two distinct 
components. The first one - virial or contact - describes pairwise interactions of grains. The second one - kinetic or 
Reynolds - is due to velocity fluctuations. 



dt 



A. Order parameter for granular fluidization: static contacts vs. fluid contacts 

At any moment of time all contacts are classified as either "fluid" or "solid" . A contact is considered "solid" if it is 
in a stuck state (Ft < ^F n ) and its duration is longer than a typical time of collision i*. The first requirements elimi- 
nates long-lasting sliding contacts, and the second requirement excludes short-term collisions pertinent for completely 
fluidized regimes. We choose a typical collision to last t* = l.li c . When either of the requirements is not fulfilled, the 
contact is assumed "fluid". 
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We define the order parameter as the ratio between space-time averaged numbers of "solid" contacts (Z a ) and all 
contacts (Z) within a sampling area p3[ . 

P {y)=Jz5fW). (31) 

This definition endures at least two limiting cases: when a granulate in a static state and when it is strongly agitated, 
i.e. c ompletely f hiidi zed. In fact, in the former quiescent state all contacts are stuck and p = 1. In the fluidized case 
(Z s ) is zero and (Z) is small but finite, therefore p — 0. 

Let us note that the order parameter p just introduced is expected to be as sensitive to the degree of fluidization as 
the stress tensor. A small rearrangement of the force network may result in strong fluctuations of either field, while 
such quantities as the solid fraction or granular temperature will remain virtually constant. It is known that granular 
aggregates exhibit rigidity phase transition near the critical volume fraction v c [ll| . When the volume fraction is near 
or above v c , the granular material has elasto-plastic rheology and below this critical value it behaves as a fluid. Near 
and above the critical volume fraction for the same boundary conditions the granular aggregate may exhibit either 
slow creeping flow, solid-like behavior, or both. The latter case, known as stick-slip phenomenon, is observed both 
experimentally and numerically 0, ITU |2S| . In an experiment under constant volume boundary conditions a series of 
stick-slip events occur without a significant change in the solid fraction, which is insensitive to global rearrangements 
relieving the accumulated stress. We expect the order parameter to be able to reflect such rearrangements and 
ultimately describe the corresponding phase transition. 



B. Stress tensor 



The full stress tensor (|30[) consists of a contact (virial) part and the Reynolds part. In turn, the contact stress 
tensor can be split into the "solid contact" component, <j s , and the "fluid contact" component, (erf), in the same 
fashion as was done with contacts themselves. Combining the "fluid contact" component with the Reynolds stress, 
we obtain the full stress tensor as a sum of two parts 

<T = a f +a s , (32) 

where summation in and Y^" is restricted to "solid" and "fluid" contacts respectively. 

The "fluid" part of the stress tensor is due to short-term collisional stresses and the Reynolds stresses, whereas 
the solid part accounts for persistent force chains. The Reynolds contribution to the stress is negligibly small in 
the vicinity of the phase transition, but comes into play when the granular aggregate is highly fluidized. In the 
system which is neither completely rigid nor completely fluidized we expect the coexistence (in time and space) of 
both phases. A particular grain may have both types of contacts at the same time thus contributing to both a* and 
a s . This picture is reminiscent of the concept of bimodal character of stress transmission in static contact network 
introduced by Radjai et al [2^.l30| . 



V. TESTBED SYSTEM: COUETTE FLOW IN A THIN GRANULAR LAYER 



A. Bifurcation diagram 

We studied the fluidization transition in a thin (50 x 10) granular layer between two "rough plates" under fixed 
pressure P and zero gravity conditions (FigurcQJ. The parameters of simulations were po = 0.5, e = 0.92, A r = 
0.2, k n = 2 • 10 5 , k s /k n — 1/3, 7„ = 16.7. We used the same material parameters of grains throughout this Section. 
The rough plates were simulated by two straight chains of large grains (twice as large as an average particle diameter) 
glued together. The layer was chosen narrow enough, so the properties of the granular layer (shear rate, shear stress 
components, order parameter, etc) were constant across the layer. Two opposite forces Fi = — F2 were applied to 
the plates along the horizontal x axis to induce shear stress in the bulk. We started with weak forces not sufficient to 
initiate a shear flow and slowly ramped them up in small increments well above the critical yield force at which the 
granular flow started. After that we ramped the shear forces down until the granular layer was jammed again. At 
every "stop" we measured all stress components, strain rate, and the order parameter, and averaged the data over the 
whole layer and over time of each step. Figure [21 shows the strain rate 7 (we drop subscripts yx at the strain rate) and 
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FIG. 1: Sketch of the granular shear flow model 



the order parameter p as functions of the shear stress a yx which is approximated by the applied force F normalized 
by the layer length L Xl for P — 40. As it is to be expected, the strain rate remains zero, and the order parameter 
is one until the shear stress reaches a certain critical value <j\ « 12.6. This value differs slightly for different runs 
because of the finite system size and absence of self-averaging in the static regime. Above the yield stress, the strain 
rate abruptly jumps to a finite value « 0.35, and the order parameter drops to ks 0.15. At larger \<r yx \ 1 the strain 
rate increases faster than linearly, and the order parameter rapidly approaches zero. The return curve corresponding 
to the diminishing of the shear stress follows roughly the same path, and then continues to another (smaller) value 
of the shear stress (era ~ 9.4). At this value the layer jams, the strain rate returns to zero, and the order parameter 
jumps back to one. 

The most striking feature of this figure is the hysteretic behavior of both the strain rate and the order parameter as 
a function of the shear stress. This hysteresis was anticipated in our order-parameter model |l8lll9j| . however now we 
are in a position to fit the model equations quantitatively. We repeated these simulations at several different values 
of the compressing pressure P. Data for different pressure values in the flow regime fall onto the same universal curve 
if one normalizes the shear stress by the pressure (see Figure |3| . Assuming that there is an (unobserved) unstable 
branch of the bifurcation diagram which merges with the stable branch at 5 m 0.26, we can make a simple analytic 
fit of this curve as 

F(p, S) = (l- p) (p 2 - 2p«p + pi cM-A{5 2 - Si)]) = (34) 

with p„ = 0.6, A = 25, (5* = 0.26 (see Figure|3 line) and use it in the polynomial expansion of the free energy density 
which enters the order parameter equation l|2U|) : 

HP) = J" F(p,S)dp. (35) 

We also measured the density and the granular temperature of the grains as we decreased the shear force. These 
measurements could only be performed in the range p < 0.55 since for larger p the partially fluidized state is unstable. 
Note that that density of grains stays almost constant in a wide range of the order parameter p > 0.1 (see Figure^). 
The granular temperature (which is defined as T = (uiUi)/2) normalized by the applied pressure P appears to be a 
unique function of the order parameter (Figure^). 

B. Relaxation dynamics of the order parameter 

We probed the relaxation dynamics of the order parameter by performing the following numerical experiment. The 
granular layer was prepared as in the previous Section. Lateral shear forces were increased adiabatically until the 
granular system reached a metastable solid (jammed) state within a hysteretic region. Then the layer was perturbed 
by applying random forces to a small randomly selected fraction of particles. The dynamics of the order parameter 
vary depending on the magnitude of the perturbation. Figure [5] shows an example of the evolution from the same 
jammed state for two different magnitudes of initial perturbation. Interestingly, the relaxation back to the jammed 
state is very fast, whereas the relaxation toward stable shear flow is much slower. We believe that it has to do with 
inertia of grains and the upper plate, so the intrinsic time scale of the order parameter relaxation is rather small, 

Unfortunately, our thin Couette flow system does not allow us to probe the local coupling of the order parameter 
since the order parameter is uniformly distributed throughout the system. In the absence of this data, this coupling 
was modeled by the linear diffusion term in (|20|l with the constant diffusion coefficient. As we will see in the following 
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FIG. 2: The strain rate (a) and the order parameter (b) vs. shear stress in a thin Couette geometry of Figure Q with 500 
particles (10x50) at P = 40. 



1.4 - 




FIG. 3: The order parameter as a function of the normalized shear stress in a thin Couette geometry. 

Section, this approximation indeed provides a good description for spatially non-uniform near-surface flow, however 
the value of the diffusion coefficient appears to be a function of local stress. 

C. Fitting the constitutive relation 

The next step is to fit the constitutive relation from MD simulations. To this end, we use the same Couette flow 
simulations, but now we analyze the "fluid stress" a a g and the "static stress" o~ s a p separately during our ramp-down 
simulations at three different values of P. Figure El shows a sample of these data for P = 30. At large 8 = F/(L X P), 
when the order parameter p is low, the total stress is dominated by the fluid component, but as the flow stops and p 
approaches unity, the fluid stress turns to zero, and the total stress is equal to the static stress. Plotting cr^ x /a yx as 
a function of the order parameter p for different P (Figure EJi) , we observe that all data collapse onto a single curve 
which is well fitted by q(p) = (1 — p) 2 ' 5 . The lines in Figure show the fit of the fluid and static stress tensors using 
iPJl.lfTr)]! with q(p) = (1 - p) 2 - 5 . 

The fluid as well as solid parts of the stress tensor are nearly symmetric, a^ x = al y s , so the ratio o~l y /a xy is described 
by the same scaling function q(p). On the other hand, the same procedure for the diagonal elements of the stress 
tensor yields a noticeably different scaling (see FigunJ7|b). Furthermore, a small but noticeable difference is evident 
between al x /a xx and Pyy/cyy More detailed analysis shows that in fact fluid parts of the diagonal components of 
the stress tensor a{. x and a^ yy are nearly identical, and the difference is due to the solid part of the normal stresses 
(see Figure |8J . This observation is consistent with the fact that the diagonal terms of the static stress tensor are 
determined by the details of the external loading. On the other hand, in a completely fluidized state the diagonal 
terms are all equal to the hydrodynamic pressure pf |3lj . In a partially fluidized regime, the diagonal terms of the 
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FIG. 4: The density (a) and the normalized granular temperature (b) versus the order parameter in a thin Couette geometry 
for three different values of pressure P. 




time 

FIG. 5: Relaxation of the order parameter towards shear flow (solid line) and jammed state (dashed line) for two different 
initial perturbations in a thin Couette system of 500 particles with P — 40, a yx — —12. 

shear stress can be expressed as 

°xx=Pf/qx(p), <7yy=Pf/q V {p)- (36) 

Both functions q x ,y{p) should approach 1 as p — > 0, but they may have different functional form to reflect the 
anisotropy of the static stress tensor. If the fluid pressure pf approaches zero in the solid state as pj oc (1 — p) a , then 
q x . y oc (1 — p 13 *'") 01 , so <J X x,yy oc (3g . In our simple Couette flow, the diagonal stress tensor components can be well 
fitted by q x (p) ~ (1 — p) 1 - 9 and q y (p) ~ (1 — p 1 - 2 ) 1 - 9 (see Figure0b). So we observe that even in a partially fluidized 
regime, the "fluid phase" component indeed behaves as a real fluid with a well-behaved "partial" pressure pf which 
is zero in a solid state at p = 1 and is becoming the full pressure in a completely fluidized state p — 0. 

Plotting the fluid shear stress versus the strain rate, we can test the validity of the Newtonian model for the stress- 
strain relation iJSJ. Figure [5] shows — o* yx vs 7 for three different pressures P = 20,30,40. At small 7 all three lines 
are close to the same straight line a^ x = 127, which indicates that the Newtonian scaling for fluid shear stress holds 
reasonably well. The deviations at large 7 are evidently caused by variations of temperature and density in the dilute 
regime. Note that in contrast, the full shear stress does not goes to zero as 7 — > (see Figure Et), so a viscosity 
coefficient conventionally defined as the ratio of the full shear stress to the strain rate diverges at the fluidization 
threshold as observed in Ref. 0. 

Combining the Newtonian law for the fluid stress-strain dependence with the order parameter scaling of the fluid 
stress tensor, we arrive at the relationship between the full stress tensor and the strain rate tensor l|19|l with p, — 
12, q x {p) = (1 - pf*, q y (p) = (1 - p 1 - 2 ) 1 - 9 , q(p) = (1 - pf*. 
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FIG. 6: Fluid and static components of the shear stress in a thin granular layer (10x50) at external pressure P — 30 as a 
function of normalized external shear stress 5 = F/(L X P): direct calculation (points), obtained from total stress using relations 
10,10 with q(p) = (1 - p) 2 - 5 (lines). 




FIG. 7: Ratios of the fluid stress components to the corresponding full stress components cr^/ua^ vs. p for three different 
pressures P: a - shear stress components, closed symbols a yx , open symbols a xy , line is a fit q(p) = (1 — p) 2 ' 5 . b - normal stress 
components, closed symbols a xx , open symbols a yy , lines are the fits q x (p) = (1 — p) 1 ' 9 , q y (p) = (1 — p 1 ) 9 . 

D. Toward dilute granular flows - granular temperature revisited 

While the above fittings have been made for the regime of a slow dense flow with v — > v c , it is tempting to generalize 
the theoretical model so it smoothly crosses over to the standard kinetic continuum theory <|Tft— (|1 Of) for p — > 0. This 
goal can be achieved by including the equation of state © and the equation for the granular temperature J2J) back 
into the theory The important difference with respect to the standard kinetic theory is that the pressure which we 
calculate with Eq.© is not the total pressure, but the partial pressure associated with the fluid part of the stress 
tensor. Of course, as p — ► 0, the static part of the stress tensor disappears, and the partial pressure becomes the total 
pressure. 

We can test the relevance of this combined approach by calculating the spatial correlation function G(y) = (1 + 
e) _1 [itd 2 pf / Ai>T — 1] with the values of fluid pressure pf, temperature T and density v calculated in our testbed 
Couette flow at different external pressures P and comparing it with the theoretical functions Gcsiy), Gcs(v), Gl{v) 
Hll[] - 113fl . Figure fTTH shows that the Carnahan-Sterling formula works very well in the dilute range v < 0.67 as 
expected. In the high density regime G(v) approaches the free volume result (|12f). and the overall dependence is in 
agreement with the global interpolation Gl(v) by Luding p4j |. 

We can carry this analysis one step further and test the kinetic theory prediction for the shear viscosity. If we scale 
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FIG. 8: Diagonal stress components (fluid, solid, and total) in a thin granular layer (10x50) at external pressure P — 30 as a 
function of external shear stress F/ L x normalized by the external pressure P. 




FIG. 9: Stress-strain rate relation for a thin granular Couette flow at three different external pressures: (a) fluid shear stress 
vs strain rate, the straight line is a constant viscosity fit a* — 12-y; inset: scaled fluid shear stress fiJ 1 cjy X as a function of 7; 
(b) full shear stress vs strain rate 



a yx by the shear viscosity calculated using the kinetic formula Q) with globally fitted correlation function Gl(v) and 
actual temperature and density values from corresponding runs, all three lines in Figure [5^, collapse onto the same 
straight line dependence fij a^ x = 7 (see Figure Et, inset). 

One more test of the kinetic theory predictions can be performed by analyzing the granular temperature as a 
function of the shear strain rate. According to Ref. ^(|, T 1 / 2 = for a plane parallel shear flow where A is a 
material constant weakly dependent on volume fraction. Our numerical data for three different external pressures 
shown in Figure ITTI are consistent with this scaling law although small deviations can be observed at small 7. 

We have not done a similar comparison for the bulk viscosity, thermal conductivity, and the energy loss, however the 
presented data strongly suggest that the constitutive relations for the fluid part of the stress tensor are well described 
by the standard granular hydrodynamics. 
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FIG. 10: Particle-particle correlation function G{y) calculated via equation of state ® using the values of fluid pressure, 
temperature, and density in a thin granular Couette flow at three different external pressures. Inset: the same data in a 
semi-log scale 




FIG. 11: Granular temperature as a function of strain rate in a thin granular Couette flow at three different external pressures. 



E. Order parameter description of partially fiuidized granular flows - Take 2 

Let us summarize the equations of the continuum theory as specified on the basis of the 2D molecular dynamics 
simulation of the thin Couette flow. 

The mass, momentum, and energy conservation conditions are expressed by Eqs. 0"©- The order parameter 
equation now has the form 

^ = D^p - (p - 1) (p 2 - 2p*p + pi exp[-A(S 2 - <5 2 )]) (37) 

with = 0.6, (5* = 0.26, A = 25. As mentioned before, for the lack of simulation data we assume linear diffusion 
coupling of the order parameter with a constant non-dimensional diffusion coefficient D. The constitutive relation 
now reads 

o a p = [Pf/q a (p) + (Pf - \f)Trj/q{p)]5 a p - p f j a p/q{p) (38) 

withg ic ( (0 ) = (l-p) 1 - 9 , g » = (l-p 1 - 2 ) 1 - 9 , q(p) = (l-pf*. 

The equation of state and expressions for viscosity, thermoconductivity and the energy dissipation have the same 
functional form as (JSJ-JTOJ, but they are now written for the "fluid" parameters pf, /i/, A/, and e. 



13 




FIG. 12: The geometry of the MD simulation of a surface-driven shear flow 
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P10V5 


5- 10 3 


50 


100 10 5 




5.0 


P10V50 


5- 10 3 


50 


100 10 50 




6.0 


P50V5 


5 ■ 10 3 


50 


100 50 5 




17.0 


P50V50 


5- 10 3 


50 


100 50 50 




25.0 


P50V5L 


1- 10 4 


50 


200 50 5 




18.5 


P50V50L 


1- 10 4 


50 


200 50 50 




26.5 


P20F10 


5- 10 3 


50 


100 20 - 


10 


20.0 


P20F20XL 


1.92 ■ 10 4 


96 


200 20 - 


20 


20.0 



TABLE I: Parameter values for the simulations of deep Couette flows for different geometries and boundary conditions. The 
first six runs were performed at constant velocity of the top plate and the last two runs were for a constant horizontal force 
applied to the plate; for all runs no = 0.3, e = 0.82, A r = 0.2, k n — 2 ■ 10 5 , k s /k n — 1/3, j n = 16.7. 



VI. SURFACE-DRIVEN SHEAR GRANULAR FLOW UNDER GRAVITY 



In this Section we apply the theoretical description which was formulated in the previous Section on the basis of 
numerical simulations of a thin Couette flow with no gravity to another model problem. We consider shear granular 
flow in a thick granular layer under gravity driven by the upper plate which is pulled in a horizontal direction, see 
Figure IT21 A similar system has been studied experimentally by Nasuno et al. Q, as well as by Tsai et al. 0. 

We simulated up to 20,000 particles in a rectangular box under a heavy plate which was moved either with a 
constant speed V or a constant force F. Periodic boundary conditions were assumed in a horizontal direction. After 
a transient, a quasi-stationary fluidization and shear flow established in the near-surface layer, while near the bottom 
grains remained in a nearly static jammed regime. Here, we show the results of several different runs with small/large 
pressure and small/large shear. The details of these runs are given in Table [I] The vertical profiles of the density, 
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FIG. 13: Density (a), horizontal velocity (b), temperature (c), and the order parameter (d) profiles in a deep granular layer 
driven by upper moving plate for four different runs from Table |H 



flow velocity, and the order parameter are shown in Figure 1131 The density of grains remain nearly constant very 
close to the maximum random packing density value, except for a narrow near-surface boundary layer. The most 
significant density variations were observed for the case P10V5Q corresponding to a light fast moving upper plate. 
The horizontal velocity decays roughly exponentially off the plate in agreement with experimental evidence 0, 0] ■ 
The vertical profiles of the order parameter demonstrate a well-defined transition from fluid state near the upper 
plate to solid state below. The width of the "interface" grows with the applied pressure P, which indicates the stress 
dependence of the diffusion coefficient for the order parameter. Surprisingly, we found that the order parameter does 
not approach 1 at large depths, but instead seems to saturate at some value slightly below 1. We believe that this 
behavior is inherent to our 2D geometry with periodic boundary conditions on side walls. The moving upper plate 
oscillates vertically and produces vibrations in the bulk of granular layer. These slowly decaying with depth vibrations 
break weak contacts between particles which are not strongly pressed against each other (e.g. lying under arches). 
We believe that in full 3D simulations with more realistic boundary conditions this effect may be less pronounced. In 
principle, it may be included in the theoretical description by proper averaging of fluctuations of the strain tensor in 
the spirit of Savage . 

By averaging velocity fluctuations and forces acting on individual particles, we calculated vertical profiles of the 
fluid and solid parts of the shear and normal stress components a a p (see Section llVjl . These profiles for Run P10V50 
are shown in Figure ITU Strong fluctuations of the horizontal component of the stress tensor a xx are mostly related to 
the static component of the tensor. According to Eqs. i|15|l . (|lt>(l . fluid and solid parts of the shear stress components 
should be related to the shear component of the full stress tensor (which in this geometry is roughly independent of y) 
via the function q(p). Figure IT31 a depicts this function as a parametric plot of a' yx {y)/<J yx vs. p(y) made using vertical 
profiles of stresses and the order parameter. As seen from this Figure, the same fit q{p) = (1 — p) 2 - 5 approximates the 
data quite well. However, unlike the zero-gravity case of the thin Couette flow, the normal stress components seem 
to be isotropic a xx — a yy and they both are well described by q y {p) (see Figure ^Jb). 

The dependence of the fluid part of the shear stress component on the strain rate (Figure I16|l shows the same 
behavior as for the thin Couette system: at small 7 the viscosity is nearly constant pf w 12, and at larger 7 shear 
thinning is observed. Interestingly, the dependence of the local Reynolds shear stress on the local strain rate is well 
described by the Bagnold scaling a yx cx j yx . Eventually, at large 7, this scaling should dominate the full stress-strain 
rate relationship. Figure IT7I compares the behavior of the viscosity coefficient p = <r yx / A / and pf = <jf x /j calculated 
along the vertical profiles of a and 7, as a function of density v and order parameter p. While the former diverges as 
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FIG. 14: Stationary profiles of the vertical (a), horizontal (b), and shear (c) stress components for run P20F10. 
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FIG. 15: Ratios of fluid and full components of the stress tensor as a function of the order parameter for four different runs at 
different speeds and pressures: a - shear stress component a yx , b - normal stress components a X m,yy Closed symbols correspond 



to a xx , open symbols correspond to a yy . Solid lines show the fits q(p) = (1 — p) J in (a), and qi(p) = (1 — p 



in (b). 
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FIG. 16: "Fluid" shear stress vs. local shear strain rate for several runs: (a) total fluid shear stress, (b) Reynolds part of the 
shear stress. 

v — > v c and p — ► 1, the latter approaches the constant value /x/ ~ 12, in agreement with results of Section Ivl 

As in the previous Section, we can extract the particle-particle correlation function by calculating G(v) = (1 + 
e) _1 [Tid 2 pf / AvT — 1] using the vertical profiles of a^ x ,T, and v. Again, we obtain a good agreement with theoretical 
predictions based on the kinetic theory for fluid component of the stress tensor (Figure EJ- 

Finally, we can compare the stationary vertical profiles of the order parameter and the horizontal velocity with 
theoretical predictions. In most of our numerical simulations we specified the velocity of the upper plate rather than 
the applied force. That allowed us to study the regimes of slow dense flows which would be unstable had we applied 
a constant shear force. The shear stress tensor component a yx in the stationary regime was indeed constant across 
the layer (see for example Figur JlU c). However, due to slippage near the moving plate the relation between the plate 
speed and the shear stress is complicated. We do not address the issue of boundary conditions here as it is a subject 
of a separate study (see for example 32] ) . Here we simply use the values which arc obtained in numerical simulations 
(the last column in Table [IJ, a s parameters in our theoretical model. 
In the stationary regime, the relevant stress tensor components are specified as follows: 

a yy = P+(H-y), (39) 

Oyx = CTxy = COUSt (40) 

where P is the external pressure applied to the upper wall. 

We also need to specify the boundary conditions for the order parameter at the top and bottom plates. This is a 
serious issue in its own right which we will address elsewhere. Here we simply impose no-flux boundary conditions 
both at the top and the bottom plate for the order parameter, d y p(Q) = d y p(H) = 0. 
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FIG. 17: Full viscosity —a yx /j (solid symbols) and "fluid" viscosity — (Ty X /j (open symbols) coefficients as functions of density 
(a) and the order parameter (b) for two runs. 




FIG. 18: Particle-particle correlation function as a function of density for several runs of the thick Couette flow simulations 



We limit ourselves with the case of slow dense flow regime, when the granular temperature plays a minor role, and 
the granular flow can be considered incompressible. This allows us to use the reduced set of equations PJl.l|37 |l .l|88 [l with 
the fixed viscosity pLf — 12. The stationary shear flow solution of the continuum equations can be found numerically 
as follows. Since the components of the full stress tensor are assumed known, we solve the time-dependent order 
parameter equation (|37l) using the pseudo-spectral method until the solution reaches a stationary state. The resulting 
solution for the order parameter is then used to obtain the velocity profile by integrating the constitutive relation Ij38(l 
from the bottom (y = 0) up. Since the grains are strongly compressed near the rough bottom plate due to gravity, 
we assume the no-slip boundary condition for the horizontal velocity at y = 0. The momentum conservation equation 
(0) is satisfied automatically. Thus obtained profiles of velocity and the order parameter were compared with our 2D 
molecular dynamics simulations. 

The results of the comparison between the velocity and order parameter profiles obtained in simulation and by 
using the continuum theory are shown in Figure ITUl for four runs P10V5 (a), P10V50 (b), P50V5 (c), P50V50 (d). 
The only fitting parameter used was the diffusion constant D in the order parameter equation, which has not been 
determined in our testbed analysis. We used D = 1 for runs P10V5 and P10V50, D = 5 for P50V5 and D = 10 for 
P50V50. From this we can conclude that the diffusion coefficient depends on the the local stress tensor, however more 
elaborate numerical experiments are needed to pinpoint this dependence more quantitatively. All other parameters 
were identical for all four cases as specified in Section IV El The vertical profiles of the order parameter and the 
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FIG. 19: Profiles of the order parameter and velocity in a thick granular layer driven at the surface by a heavy moving plate 
for Run P10V5 (a), P10V50 (b), P50V5 (c), P50V50 (d). Lines show the theoretical results obtained from the continuum 
model (l37L 11381 . empty symbols indicate numerical data. Insets show the velocity profiles in the logarithmic scale. 



horizontal velocities are reasonably well described by the theory. However, for low pressure runs P10V5 and P10V50, 
the horizontal velocity profiles deviate from the numerical data presumably because the viscosity coefficient is no 
longer a constant in a dilute region near the top plate. 



VII. CONCLUSIONS 



In this paper we performed a series of numerical simulations of 2D wall-driven granular flows. These simulations 
were designed with a specific goa l, to quantify the continuum theory of partially fluidized granular flows which was 
introduced ad hoc earlier [18|, [13 • We defined the order parameter as a ratio of the number of static contacts to 
the total coordination number averaged over a small mesoscopic volume. Using simulations of a thin Couette flow 
between two rough plates, we determined the free energy density for the order parameter. Simulations confirmed that 
the ratio of the shear to the normal stress in the bulk of the granular flow can parametrize the stationary states of the 
order parameter equation. The same simulations allowed us to determine the detailed structure of the constitutive 
relation. We split the total stress tensor into the fluid and solid components, in which the former is comprised of the 
Reynolds stresses and the stresses transmitted through short-term collisions, while the latter is formed by the force 
chains through persistent contacts. The ratio of fluid and solid stress components is indeed determined by the order 
parameter through scaling functions q{p), q x . v (p). Remarkably, the fluid component of the stress tensor is a linear 
function of the strain rate 7 in the slow dense flow regime. This justifies the Newtonian scaling of the stress-strain 
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relationship adopted in the theory. 

Using the calibrated theory, we studied the flow structure of a thick surface driven Couette granular flow under 
gravity. We found the theoretical predictions to be in a good quantitative agreement with simulations. 

The evidence presented here suggests an intriguing interpretation for the order parameter description of dense and 
slow granular flows. The granular material under shear stress appears to be similar to a multi-phase system with the 
fluid phase "immersed" into the solid phase. The fluid phase behaves as a simple Newtonian fluid for small shear 
rates when the density is almost constant, but exhibit shear thinning at larger shear rates when the density begins to 
drop. We observed that the Reynolds part of the fluid shear stress obeys the Bagnold scaling a^ x ~ j 2 . We anticipate 
that for very large shear rates when the Reynolds stress becomes dominant the overall stress tensor should exhibit 
Bagnold scaling locally. 

While this theory is primarily intended for dense and slow granular flows, we have shown that it can be combined 
with existing models of rapid granular flows based on the kinetic theory of granular gases. This requires to drop the 
assumption of incompressibility and include the equation for the granular temperature. Our simulations showed that 
the kinetic theory works well for the fluid part of the stress tensor in the whole range of densities from dilute regime 
to the critical random close packing density. 

Many issues still remain open. The spatially non-uniform dynamics of the order parameter requires a more detailed 
study. We found that the diffusion constant postulated in Eq. i|37|) appears to be a function of the normal shear stress 
as well as the local strain rate, however, we do not have sufficient numerical data to provide a quantitative description 
of this dependence. It would be of interest to analyze the propagation of a fluidization front in a granular layer 
prepared in a meta-stable static regime. Such simulations could provide an insight into the mechanisms of the local 
coupling of the order parameter. 

The molecular dynamics algorithm employed is based on a number of approximations. These approximations, 
however well tested and widely accepted [H El ^3 i directly affect the results of our fitting the continuum model. For 
example, if one replaces the Hookian model of particle interaction with a Hertzian one, an appreciable difference in 
the structure of the order parameter may be observed. More numerical work is needed to quantify the relationships 
between the microscopic parameters of the system (nature of collisions, restitution coefficient, friction, etc) and the 
parameters of the continuum model. 

Finally, our simulations were limited by 2D systems, and of course the resulting continuum theory can only be 
directly applicable to 2D systems. While we anticipate that the structure of the model to remain in 3D systems, the 
specific form of the fitting functions should change. This future work will allow us to perform a comparison of the 3D 
model not only with numerical simulations, but also with experimental data. 

The authors are indebted to B. Beringer, P. Cvitanovic, J. Gollub, T. Halsey, and J. Vinals for useful discussions, 
and to J. C. Tsai for sharing his unpublished experimental data. This work was supported by the Office of the Basic 
Energy Sciences at the US Department of Energy, grants W-31-109-ENG-38, and DE-FG03-95ER14516. Simulation 
were performed at the National Energy Research Scientific Computing Center. 
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